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f-H I A topological approach and understanding to the detection of unstable periodic orbits based 

'~{ . on a recently proposed method (PRL 78, 4733 (1997)) is developed. This approach provides a 

classification of the set of transformations necessary for the finding of the orbits. Applications 
, to the Ikeda and Ifenon map are performed, allowing a study of the distributions of Lyapunov 

exponents for high periods. Particularly the properties of the least unstable orbits up to period 36 
are investigated and discussed. 
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^ I. INTRODUCTION 

Unstable periodic orbits (UPOs) represent a skeleton of complex chaotic systems and allow the calculation of many 
characteristic quantities of the underlying dynamics like Lyapunov exponents, fractal dimensions and entropies of the 
attractors . For dissipative systems expansions in terms of periodic orbits are well established in the literature 
and demonstrate the relevance of the cycles for the understanding of chaotic dynamics. Both low-dimensional 
\^ model systems such as discrete maps as well as experimental time series ||7|-p^ have been studied. Furthermore 
I the series expansion of semiclassical properties of classically chaotic Hamiltonian systems with respect to the length 
• and stability coefficients of the periodic orbits is a fruitful and frequently applied technique: it allows the investigation 
' of the energy level density as well as other quantum properties |ll[|. Much has been added to the importance of the 
^ UPOs by using them to control chaotic dynamical systems ( see |12|] and references therein ). 
• '-H ' More recently, cycle expansion techniques have been invented and shown to converge well, especially when the 
^ , symbolic dynamics is well understood [^,^. Series expansions over periodic orbits used for calculating dynamical 
' averages are typically ordered according to the orbit length p ||l|,p^-p^ . Drawbacks of these expansions are the large 
. ^ , number of orbits (increasing exponentially with p) , the required completeness of the set of cycles for given period and 
' the slow convergence properties ||l^,|l5|. A promising proposal was made [ p^|JT^] saying that series expansions could 
H [ converge better if they are truncated according to the stability of the orbits ]18[ . What is more, stability ordering 
. . . 1 does not rely on the knowledge of the systems symbolic dynamics, which is unknown for generic dynamical systems. 

Chaotic dynamics is intrinsic for many physical systems, which is why periodic orbit theory is not restricted to 
specific areas of physics. Much effort has been spent on developing efficient techniques for calculating UPOs of a 
given dynamical system. What makes this quest so difficult is the exponential proliferation of the number of cycles 
with increasing period and their increasing instability (in fully chaotic systems). Specifically the modified Newton- 
Raphson method is an all-purpose method, which does not require a special form of the underlying equations of motion. 
However, it needs a good initial guess for the starting points. Therefore this algorithm rapidly becomes expensive and 
is limited to relatively short periods and low-dimensional systems. A variety of other methods have been developed 
which focus either on time series analysis or on finding UPOs for given equations of motions [pP Jl^ , p^| -pl| . For a 
special class of systems, a numerical technique for calculating arbitrarily long UPOs to any desired accuracy was 
introduced in |2^] for the Henon map |^ and later applied to a few other dynamical systems [ p^p^p^ . This method 
allows the systematic computation of all UPOs of any given order, each given by a unique binary symbol sequence. 
Recently, we proposed an alternative method (in the following referred to as SD method) for the calculation of 



UPOs |25|. The basic idea is to transform the fully chaotic system to a new dynamical system with the periodic orbits 
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keeping their positions but changing their stabihty properties: For a particular transformed system a certain fraction 
of the periodic cycles becomes stable and can be found by simply iterating the transformed system. This fraction 
depends on a tuning parameter which represents an upper stability cutoff for the fixed points (FPs) to be detected. 
In the following we will use the term 'stabilisation' and 'stabilised fixed points' for the process and fact respectively 
that the unstable FPs of the original chaotic dynamical system have become stable in the corresponding transformed 
system. The reader is kindly asked to distinguish this use of the term 'stabilisation' from the one used in control 
theory of chaos via unstable periodic orbits. According to the above the SD method allows the systematic calculation 
of the least unstable periodic orbits of any given order p [p6| . The latter possibility meets the requirements of the 
series expansions using stability ordering, since they allow to derive properties of a physical system by exploring only 
the least unstable period orbits up to a given stability cutoff. The technique of the SD method is highly flexible 
and can be applied in a straightforward manner to a great variety of discrete dynamical systems of any dimension. 
Continuous dynamical systems can thereby be treated using properly chosen Poincare surfaces of section. Therefore, 
it might open the way for employing the proposal of p^ , p^ for a great variety of dynamical systems. 

The present work has a twofold aim. The first goal is to enhance our understanding of and give insights into the SD 
method. To that end we provide a complete classification of the involved transformations. As a result we gain both 
topological as well as geometrical understanding and interpretation of the transformations. Corresponding invariant 
structures are thereby revealed and the FPs can be classified similarly to the stabilising transformations. This opens 
the future perspective to selectively detect UPOs not 'only' with respect to their stability but also with respect to 
certain desired geometrical properties. We will thereby learn how simple global operations on the dynamical system 
change the stability properties of fixed points. The second aim of this work is to elucidate and extend the work in ref. 
psf . To this end we provide extensive results of applications of our approach to the Ikeda and Henon maps. 
For the Ikeda map, we calculate the complete sets of FPs for periods up to p = 15. The number of orbits is large 
enough to investigate the distribution of Lyapunov exponents. This distribution is compared with the corresponding 
one of the Henon map according to . Next the algorithm of to determine the least unstable orbits is slightly 
modified and applied to the Ikeda and Henon map. The ten most stable orbits of a given length up to period p = 36 
are studied for both maps. Their Lyapunov exponents vary in a characteristic way as a function of the period. 

In detail we proceed as follows: In Section 2.1 we briefly review the SD method for finding UPOs as described 
in The topological/geomctrical classification and geometrical extension is presented and discussed in Sec. 2.2. 
In Sec. 3.1 we show the results of complete sets of orbits and we present and analyse the distribution of Lyapunov 
exponents in Sec. 3.2. The extended method to stabilise the least unstable cycles is given in Sec. 3.3 together with 
the corresponding results for the Henon and Ikeda map. Sec. 4 concludes with a summary. 



II. THEORY OF THE STABILISATION OF FIXED POINTS 



A. Brief Review of the underlying Method 

In order to be self-contained and for our further theoretical investigation we recall in the following the key ideas 
of the method developed in ref. to detect UPOs. Consider a discrete chaotic dynamical system U : r^+i ~ /(fi) 
in n dimensions. The FPs of the pth iterate J^^' are points of the UPOs of period p. To find the FPs of U the 
following strategy is employed: A set of transformations is specified which transforms U into new dynamical systems 
{Ska} with the FPs keeping their original locations in space. The set {Ska} is chosen such that for each unstable FP 
au of U there exists a specific transformed system Sk'a' of the set {Ska} for which this FP has become dissipatively 
stable (i.e. Re(Ai) < for continuous systems, |Ai| < 1 for discrete systems, with eigenvalues Aj) and can therefore 
be detected by simply iterating properly chosen starting points of the transformed system Sk'a' ■ For each Sfa' a 
different set of FPs of U is stabilised. Let us denote by /min the minimal set of pairs (ka) with the property that there 
exists for each unstable FP a„ at least one pair {k'a') G /min for which Sk'a' transforms a„ into a stable FP as- This 
set holds for arbitrary period p. The search for the FPs of U is then straightforward: A starting point chosen in the 
global neighbourhood (see below) of the FP iterated with the transformed dynamical system Sk'a', {k'cr') G /min 
converges, due to the stability of a^, to the position of as which is equal to that of a„. 

Propagating a set {vi} of starting points and using all {ka) G /,nin we end up with a set of FP of U whose 
completeness can be ensured by successively enlarging the set {r,}. Let us now specify the systems Ska'- 

Ska ■■ Tj+l + Aka[f{ri) - Yj] (1) 

Ska are linear transformations of the original dynamical law U . Aka are invertible constant n x n matrices. The 
definition of Ska and U clearly shows that their FPs are one to one and at the same positions. Eqn. (|l|) implies the 
following relation for the stability matrices Tu and Ts^^ of U and Ska , respectively: 
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Ts,„ =l + Ak.(Tu-l) 



(2) 



In ref. it was shown that Aka can be cast in the form = A ■ Cka with < A < 1. The set of matrices {Cka} 
contains aU orthogonal matrices with only one non- vanishing entry ±1 per row or column, i.e. they represent a group 
of special reflections and permutations. 

(T = ± indicates the sign of the matrix determinant, and k is an additional label to uniquely specify the matrices. 

In two dimensions we have Co+ = — C2+ — {qi), Cq- = — C2- = {~q1), Ci+ = — C3+ = (-"lo) ^1- = 

— C3_ = ~g^^ ■ The matrices {Cka\k = 0...3; a = ±} form a group with {Ck+\k = 0, 3} (matrices with positive 

determinant) being a subgroup of order 4. Table ^ is the corresponding multiplication table. Obviously, the product 
of two matrices is 

Ck"a" = Cka ■ Ck'a' with k" = k' + a'k mod 4 and a" = aa' (3) 

The notation introduced above is different from the one used in ref. [ p5| and will reveal its meaning in the course 
of the classification of all possible FPs (occurring in the original and transformed system), as provided later on. 
The arithmetics with respect to the first index {k, k' , k", ..) has always to be taken modulo 4. We remark that the 
minimal set /min is significantly smaller than the set of pairs (fccr) belonging to the matrices {Ska}- Given a certain 
unstable FP the above choice of linear transformations represented by the set {Cka} of matrices allows to find a 
particular Ck'a' such that the real parts of the eigenvalues of Ck'a' {Tu — 1) are negative at the position of the FP. As 
a consequence (see Eqn. (||)), if A is chosen sufficiently small, the magnitude of the real parts of the eigenvalues of the 
FP as in the transformed system are smaller than one and we therefore encounter a stable FP which can be detected 
as described above. The critical value of A, which just suffices to make the FP stable, can immediately be read off 
from the quadratic equations relating the stability coefficients of the original system U and those of the transformed 
system Ska The above procedure can easily be extended to higher iterates /'^'(r) of U (by replacing / with /'••^^ 
in Eq. (^) allowing us to determine all order p cycles of U. 

The advantage of the SD method is clearly its global character in the sense that even points far from the linear 
neighbourhood of a FP are attracted close to the FP after a finite number of iterations of the transformed dynamical 
law. The basin of attraction of a single stabilised FP is a simply connected area in phase space. The typical number of 
starting points needed to obtain the UPOs of a given length on the attractor is only slightly more than the expected 
number of cycle points themselves. 

The parameter A is a key quantity here. It is related to the stability of the desired cycle in the transformed system. 
With increasing period of the cycles, A has to be reduced to achieve stabilisation of all FPs. One may, however, also be 
interested in the most stable periodic orbits of a given period p p6| which is one of the key issues of the present work. 
In this context A operates as a filter allowing the selective stabilisation of only those UPOs which possess Lyapunov 
exponents smaller than a critical value. Therefore, starting the search for UPOs within a certain period p with a value 
A = O(10~^) and gradually lowering A we obtain the sequence of all unstable orbits of order p sorted with increasing 
values of their Lyapunov exponents. In |^ it was shown that for the specific choice (ka) = (0+) G /min the relation 
between A and the stability coefficients of the FPs of the original system [/ is a strict monotonous one. Transformed 
dynamical systems Ska belonging to other pairs (k'a') ^ (0+) do not obey such a strict behaviour but show a rough 
ordering of the sequence of stability eigenvalues of the FPs of U stabilised in the course of decreasing values for A. 

B. Geometrical Interpretation of the Stabilising Transformations Sku 

1. Classification Scheme 

The stability properties of the FPs of the dynamical systems Ska have been investigated in refs. |^ exclusively in 
the context of their relation to the stability coefficients of the original system U. To gain a deeper insight into the 
geometrical meaning and the interpretation of the transformations Ska which turn unstable FP into stable ones, one 
has to go beyond the pure consideration of their eigenvalues at the positions of the FPs. 

In the following we develop a geometrical approach allowing us to classify the FPs which are stabilised by different 
matrices Cka- We will hereby focus on systems with two degrees of freedom. Expectedly there should, however, 
be no major obstacles with respect to the generalisation to arbitrary dimensions. When dealing with the stability 
transformations the natural problem arises: Restricting ourselves to the set of orthogonal Cfccr-matrices with exactly 
one non-vanishing entry (±1) in each row and column and to the linearised dynamics around a FP, what can we say 
about the action of the matrices Cka (see Eqn. (|l|)) on this simple dynamical system? To approach this problem 
consider the following set of equations: 
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x = F(x) , i?(x) = /»(x)-x , F(x) = (Fi(x),i^2(x))^ 



(4) 



which describes a vector field around the FP located at x/, where F{xf) — (the superscript denotes the trans- 
posed). In the following sections we generally focus on the discussion of the system F{xf), unless noted differently. 
Now we apply the transformation 

X = FkA^) = Cka ■ F(x), FkA^) = (Ffc,4(x), Ffe,,2(x))^ (5) 

It is important to note that the dynamical system Ska in Eqn (0) represents a discretisation of the continuous system 
. Multiplication with Cka intermingles the x~ and y-coordinates of F(x), which in general changes the eigenvalues 
and eigenvectors of the corresponding stability matrix. One direction of the above problem is: Are there any points 
in the neighbourhood of x/ where this change of the dynamics is controllable? In fact let us consider the manifolds 
Zi, Z-2. defined by |3 

= {x|F,(x)=0}, i = l,2 (6) 

In the linear neighbourhood of x/ these sets clearly define straight lines. In the more general case of a nonlinear 
system they are implicitly defined continuous curves in an area of phase space where the Jacobian of the map is 
regular. Their intersection is Z\V\Z2 = {x/}. With Cfcg. acting on Zi, Z2, they either stay the same if Cfccr does not 
interchange the coordinates, or they are transferred one onto the other if Ck„ does interchange the coordinates. In this 
sense, the manifold Z = Zi U Z2 is invariant with respect to application of the set of matrices Ckin i-c Cka-{Z) = Z 
for all (ka). 

In the following we derive a classification scheme for the linearised dynamics around a FP whose validity is however, 
due to the global character of our approach, not limited to the linear regime. In order to elucidate the action of the 
stabilising transformations and to distinguish between FPs with different stability properties, let us introduce two 
different ways of classifying the matrices of a two-dimensional dynamical system, each providing its own insights. 
The classifications are such that they reflect certain geometrical features of the flow around the FP. These features 
are the different invariant sets Z on the one hand and additional geometrical properties of the matrices which have a 
particular Z in common on the other hand. 

The flrst classiflcation introduces classes consisting of matrices which have the manifold Z ^ Zi U Z2 itl common. 
They are labelled C(0,nin, ^max), where (j)min and (j^max sue the azimuthal angles of the manifolds Zi and Z2, respectively, 
being sorted with increasing order, (fimin and (/)max are related to the stability matrix of the FP in the following way: 
The linearised dynamics of eqn. (0) in the neighbourhood of a FP reads x = B • x, where B = (aij) i<i<2 is the 

l<j<2 

stability matrix of F{x.) at the FP and x is the displacement with respect to the FP. It can be shown that (f>rain and 

max 

are given by 

((/'min,(/'max) = (niin{(?!)i},max{(?!)i}) with 0i = arctan ( - — ) , i = l,2 (7) 
\ I I ) \ ai2 } 

If a stability matrix B belongs to the class C(0i„inj 0max) also its products {Cfco- • B} belong to this class. 

For the later on discussion we introduce here three sets of FPs each of which is an infinite unification of classes 

^(0mini 0max) ■ 

^ — {^(0mini 0max 

)l < 

mini 0max 

<7r/2} 

)| < 0min < 7I'/2 < 0niax < Tt} (8) 
= {C(0min, 0max)| 7r/2 < ^^jn, ^max < Tt} 

A further classification of the matrices within each class C(0min, (?!>max) is needed for a more detailed identification of 
the geometrical properties of the flux around a particular FP. To this aim, we assign a label (It) to each stability 
matrix B;,- of FPs with the following meaning: r — ±1 gives the sign of (detB;,-). To illustrate the meaning 
of /, we write i^(x) = (rcos?/', r sin?/')"^ in polar coordinates and consider the azimuthal angle V'min of -F(x) for 
X — (cos (/)min, sin (/imin)"^. By Construction, 0,nin can take the values rmr/i^ m = 0, ...,3. Now we define the index 
I = m + T ~ 1 mod 4. Encircling the FP on a unit circle x = (cos (f>, sin 4>)'^ , the normalised flux F(x)/|F(x) | describes 
a circle in the local coordinate system centred in x, too. r gives the orientation of this circular rotation of the flux 
(r = +1: anticlockwise, r = —1: clockwise) whereas / is directly related to the phase of the flux at (jimin- This 
naturally introduces a subset Cir{(t>min, 4'miiyi) of the class C(0i„in, ^max): those members of C(0,„i„, </)i„ax) belong to 
C(r(0min, '/'max) which posscss the indiccs (It), i.e. the sign r for the rotation of the flux and the phase I of the flux at 



4 



0min- We can now allow 0min, '/'max to vary while keeping the indices (It) of the matrices of this set fixed. For each 
(It) we thereby form a class Aw = {Cir{4>min, </'max)|0 < (/jmin, <^max < 27r}. For fixed 0,nin, 0max, multiphcation by the 
matrices Cka transfers one complete set Air into another set Ai'r'- The corresponding transitions are given in table 
||. The asterisks in the first three columns of table || indicate the sets C* and Air to which stability matrices of a 
FP of a two dimensional chaotic system can belong. This is of relevance when asking for the possible sets of matrices 
Cka which stabilise all FPs of a dynamical system (see below). As one can read off table |l|, the matrix Cka necessary 
to transfer stability matrices of a class Aw into matrices in another class Apr' is given by 

Cka '■ Air — > Ai'r' with I' = k + (T mod 4 and t' — a jr (9) 

It is important to note that table || holds not only for the unified sets Aw but particularly also for the individual 
subsets C/T-((/)inin, (^niax)- It IS a remarkable result that the multiplication table of the matrices Cka, i-e. table 0, 
considered as a transition table, has the same entries as the transition table ||. The multiplication law eqn. (|^) 
becomes a law for the transition of C^a to Cy a" via Cy a' with the corresponding indices k' — k" ~ aa" k and 
a' = aa" . This suggests that the matrices in a class Aw are in a way similar to the matrices Cka with k — I and 
cr = r as far as the rules for multiplication with the matrices Cka are concerned. 

The class C/r (0min, '/'max) coutaius still an infinite number of (stability) matrices. However, to gain relevant infor- 
mation on the stability properties of FPs it suffices, as we shall see in the following, to know to which of the sets 
{C^\i — 1...3} and additionally Aw the stability matrix of the FP corresponds. 



2. Properties of the Angular Functions ip[cj>) of the Flux and Examples 

Let us consider an arbitrary but fixed stability matrix B;^ G C((/)inin, (^max), ^w G Aw, and form the matrices 
Cka • B;,-, fc = 0, ...,3;(7 = ±1. Each of these matrices is element i.e representative of a different class Aiw' and we 
therefore have 

B/v = Cfcg. • B;^. with I' = k + al mod A and r'^a/r (10) 

In the following we will call the set {B;t-|/ = 0, ...,3;r = ±1} the family of the matrix B;,-. A central issue is to 
analyse which members of the family of a FP are stable. Once the information on the stability of the members of 
a family becomes available this represents an important step towards the selective use of eqn. (|l|), i.e. the selective 
detection of FPs. 

In the following we study the oricntational properties of the fiux for a family of stability matrices with fixed but 
arbitrary (0min, 0max)- For the matrices B;,- = {i'^ip)i<i j<2 ^ Aw we introduce the angular functions ipwi4') of the 
flux X = B/t- • X at a point x = (cos cj), sin (f>)'^ 

ipw{(l>) = arctan + 7^ f | - sign(B;^ • x)2) (11) 

sin ( 



(Bz, -x)! V2 

/ 6,7^ cos (h + bl'j"^ sin \ / 3 . /, (w) , , (w) 

= arctan ^2 Z + 7^ - - sign bX^' cos 6 + br > 

\b^lpcos<f> + b^l^hm^J V2 ^ V 21 ^ ■ 

il^w{4>) is the azimuthal angle of the fiux (given by F;7-(x)) with respect to a local polar coordinate system centred at 
the displacement x = (cos 0, sin (j))^ . This is illustrated in Fig. ^ Let us derive relevant properties of the angular flux 
functions ipir{(j))- i!ir{4>) is a continuous function of (p with 4'It{4') £ [0, 27r]. Due to symmetry reasons it is sufficient 
to consider the range < cj) < n. Furthermore 'ipir{(f>) is defined mod 27r and ipw{0) — '0ir(7r) = tt. Its derivative reads 

/' det(B;^) 
LwW 

with 



Lir{<P) = |x| = y{b''^^hos(f> + b^^2^sin(l>)^ + (6^\'''cos(/) + 6^2''^sin(/))2 (13) 

being the "velocity" of the fiux in x. Since multiplication with any Cka affects complete rows of B;,- by interchanging 
them or inverting their signs, Lf^{(f>) is invariant, i.e. the same for any of the resulting matrices B/t-. Two functions 
i^wW and il'vrii') therefore differ only by a shift: ipw{4') ~ '^i't{4') = [l — I' mod 4) • / can be interpreted (see also 
subsection 2.2.1) as the phase of the linearised flux around the FP: 
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V'ir(0mi„)-a+l-r).| (14) 

Accordingly we have det(Cfco- • B) — det(Cfccr) • det(B), with det(Cfe_|_) — +1 and det(Cfc_) — —1 and therefore 

V^,V(0) = -^PLM (15) 

which yields 

llJlr{(f>) = 2 ■ i>ir{(l)min) - Ipl-rift)) (16) 



T = ±1 indicates the sign of det(B;T-) and determines whether iI)it{(I>) is rising or falling, according to eqn. ([121). 

In the remaining part of this subsection we provide generic examples and illustrations of the linearised dynamics 
around a FP and the corresponding angular functions of the flux for the different members of a family {B;r}. 

The stability matrices B/,- in Figs. ||a) and Fig. |^a) are members of a family of the matrix Mi = (^_\o ^^)' 

corresponding to the class C(0.24, 1.10). Figs. ^) and |5^) are obtained in the same way, showing the family of 
M2 = ( 10 1;"') , which is in the class C(0.24, 2.03). Each subfigure shows the linear neighbourhood of the corresponding 
FP. {x,y) are the coordinates of the displacement of a trajectory with respect to the FP. The manifolds Zi, Z2 are 
displayed as long dashed lines with the direction of the flux on the Zi, Z2-lines being indicated with open arrows. 
The thick lines with flUed arrows show the position and stability properties of the eigenvectors (For saddle points, Ue, 
Vc are indicated by outward and inward looking arrows corresponding to the expanding and contracting manifolds, 
respectively. For sinks and sources the corresponding eigenvectors are labelled We,i, Ve,2 and Wc,i, ^''c,2, respectively). 
Additionally, some trajectories have been included to visualise the direction of the flow around the FP. In Figs. ||b) - 
||b) the corresponding angular functions V';t(</>) are plotted. It is important to note that these particular functions - 
as well as the phase portraits given in the corresponding flgures - are merely examples to demonstrate the qualitative 
variation of the dynamics generated by members of a family of stability matrices. There exist in general other families 
in the same class C(0„iin, i'mayi) H Ait whose phase portraits and functions ^-)ir{4>) look different from those displayed 
in the figures. However, the functions il)ir{(j)) of all these matrices in C(0inin, '/'max) H Ait have the values of il)ir{(j)) in 
4> — 0min and (j) — (/)max and the sign of ili'{(j))iT in common. The actual ipiT^fj}) can vary in between according to the 
definition (111]) . Nevertheless the above information about the stability matrices is in fact sufficient to determine their 
stability properties. 

Now, looking at the examples in Figs. ||a) - ||a) except the cases of spiral points Bi+ and B3+, we are already in the 
position to suggest a criterion for the stability and the approximate position of the eigenvectors of the FP. Consider 
the direction of the flux in each of the sectors which are determined by the manifolds Zi, Z2 and the coordinate axes. 
The manifolds Zi, Z2 are by deflnition the lines for which the flrst, respectively second component of the flux vanishes. 
Since the flux -F(x) is a continuous vector function of the angle (f> there have to be certain values 0e for which the 
flux F{x.) is coUinear with the position vector x = (cos </>, sin c/))-^, i.e. F{x.) = ±c-x, c > 0. These angles (/)e obviously 
are the polar angles of the eigenvectors of the corresponding stability matrix. The intervals where these angles 0e are 
located are bounded by Zi and Z2. The manifolds Zi and Z2 represent the boundaries of sectors in which the angle 
of the flux varies by it/2. For reasons of continuity of iI'It{4') { see Figs. |^)-|^)) one or even two eigenvectors are 
located within these sectors, which are shaded grey in the flgures. In specific cases parts of the sectors defined by Zi 
and Z2 are excluded for the localisation of the eigenvectors since coUinearity of x and F(x) is not possible. In these 
cases the outer coordinate system (x, y) represents the boundary of the sectors for the occurrence of coUinearity. 

The fact that 4>e is the polar angle of an eigenvector of the stability matrix A is equivalent to 

I (± \ _ \ 4'e ■ Re(A) > 0, unstable eigenvector , 
^i^(0ej - : Re(A)<0, stable eigenvector ^ '^ 

This means that a crossing of ^ir{(j)) with 

Xn{(t)) = + n-7r, n = 0, 1 (18) 
indicates an unstable or stable eigenvector for the corresponding value of </>, respectively. 



3. Stability Properties of the Classes and the Alternative Sets of Stabilisation Transformations 

We will show in the following that for fixed (pmin and 0max two complete classes Ai't' and Av^t" are related to 
stable FPs, i.e. the matrices in these classes have eigenvalues Ai, A2 for which Re(Ai), Re(A2) < holds. To these 
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two classes correspond two matrices Ck'a' and Ck"cr" which transfer an original stability matrix in a class Air into 
the desired stable matrices in Ai'r' and Ai"r"- This corresponds to a transformation of the original FP into the 
desired stable ones in the transformed dynamical systems S'fev and Sk"cr"- The two matrices Ck'a' and Ck"a" which 
accomplish this transformation can be read off the transition table ^ immediately. It will in particular become evident 
that a minimal set of three Cka matrices is sufficient to stabilise any FP of a two dimensional fully chaotic system. 
What is more, this classification proves to be useful not only for saddle points but also for systems with repellors and 
spiral points, which can be transformed to stable FPs (sinks) via a certain matrix Cka (see below). 

We now discuss the properties of all possible stability matrices according to their ipi^ {(j>)-dia,giams. To do this it 
is sufficient to consider the assignment of the matrices to the classes C^, and as introduced in (||) and to the 
classes Air- In the following we concentrate on the classes and since the argumentation for the class is 
analogous to that for C^. For the same reason we restricted the examples in Figs. to the classes and C^. In 
Figs. Hc)-||c) we show the areas in the ■0;i-((?!')-diagrams where a crossing of ipiT{(t>) and Xo(0) or Xi(<^) niay occur as 
gray shaded boxes. They are derived by simple application of continuity arguments concerning tpir{4>)- The labels It 
in the boxes are the labels of the corresponding class Aw of matrices whose real eigenvectors have azimuthal angles in 
this range of (f). Two labels given in brackets indicate the possibility of either two real eigenvalues with eigenvectors 
in this range (sink or source) or complex eigenvalues without real eigenvectors of the corresponding matrix (spiral 
points). These two possibilities cannot be distinguished within our classification scheme of matrices. However, this 
fact does not affect the final selection of the minimal stabilising set of matrices Cka- 

We begin our discussion with stability matrices B with negative determinants: 

• Matrices B e {^/-|^ — 0...3} are stability matrices of saddle points. The stability properties of these FPs 
are easy to determine from the corresponding ipiT{(j>)- diagrams in Figs, [sja), ^). It is obvious from the 
monotonicity and continuity of the '(/'ir ('^)-curve that they intersect the "unstable" and the "stable" lines Xo(0)i 
Xi ((/)), respectively, exactly once. The sectors where the corresponding eigenvectors are localised are shaded 
grey in the corresponding phase diagrams Figs. |^c), ^c). 

FPs with positive Jacobian B are a bit more subtle matter. They belong to the classes {-4;+|^ = 0...3} and include 
sinks and sources as well as spiral points for which the real part of the stability eigenvalues Re(Ai) Re(A2) possess the 
same signs. Spiral FPs possess stability matrices B with eigenvalues Im(Ai), Im(A2) ^ 0. This implies det(B) > 0. 
As pointed out earlier the classification of the stability matrices with C(0min, ^max) and Aw does neither specify 
the matrix nor its eigenvalues completely. In particular for stability matrices whose family contains spiral points an 
additional criterion is needed to analyse their stability since the corresponding ■!/';T(</')-functions do not cross the lines 
Xo(0) or Xi (</')• For this analysis we suggest the following criterion (whose completeness we could not prove yet). 

Consider the angles (j)t, t = 1,2, for which ■il;[^{(j)t) = 1, i.e. 4'It{4') is tangential to x«(0)j and define the distance to 
both lines Xn(0) 

dn = min |'(/';T(</>t) - </>t - f^7r| , n = 0, 1 (19) 

Now the stability of the spiral point is suggested by 

do < di => Re(Ae) > 0, unstable spiral point (20) 
do > c?i => Re(Ae) < 0, stable spiral point (21) 

This implies that for spiral FPs the line Xn{4') which is closest to i^iri't't) determines the stability and this can be 
seen as a generalisation of the criterion of the crossing with Xn(0) in the case of real eigenvalues. 

It is an interesting property that the stability of matrices B;+ of a family {B;t-|Z = 0...3,t — ±} changes with 
1^1 + 2 mod 4, i.e. if Bi+ is stable, then B3+ is unstable, if B2+ is stable, then B4+ is unstable and vice versa. 
If the eigenvalues are real the eigenvectors are not affected by the shift (only their stability changes). This property 
is obvious taking into account that the shift in I corresponds to a shift of ipirifp) by tt, with tt being also the shift 
between the lines Xo (</'): Xi(0)- Physically, the shift by tt is for any stability matrix equivalent to a time reversal 
(which corresponds to a reversal of the flux F(x) — > —F{x). 

We now come back to the discussion of the stability properties of matrices B possessing positive determinant: 

• Matrices B e — 0...3} describe sinks and sources and spiral points. We first address the sinks and 
sources and second the spiral points. The corresponding curve tpiT{4>) crosses one of the lines Xo('/') (sink) or 
Xi(0) (source) twice. It is obvious but nevertheless important to note that there cannot be more than two 
crossings, which can also formally be seen from eq. (p^). Having three or more cuts with Xn('?^), '^ = 0, 1 implies 
that 'ipir{<p) has at least two turning points in any interval [(pd, 4'u] with (j)d — 4>u = 7r/2. This means that Lf^{(f)) 
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has two or more extrema in [4>d,4'ii]- Differentiating Lf^[(j)) (see eqn. (p^) yields (for tfie sake of simplicity we 
omit the superscripts (It) in the entries ibf,^^] of B) 

^ ^I'^^f'^ = -{bj-. + 6?,)2cos0sin0 + (bL + 612)2 sin cos + 
d (p 

2(cos2 (/) - sin^ <p){b2ib22 + 611612) 
= (^'22 + &?2-&2i-&n)sin(20) + 

2(cos2 (j) - sin^ 4>){b2ib22 + 611612) cos(2<?!>) (22) 



d H2 + bh - Hi -bli 

which has exactly one solution for (j) in any interval [(f)d, 0u] with (pd — <pu = ''■/2, i.e. there is only one turning 
point in this interval. 

— The class C^: As can be read off directly from the diagrams Figs, ^jb) and c), matrices in fl A2+ and 

n ^2— &re sinks and sources, respectively. One of the eigenvectors is in the interval [^minj^max 

1, the 

other one is located in [7t/2, tt]. Matrices of the class nAi+ and n^3+ are either sinks and sources or 
spiral points. For real eigenvalues matrices of n^i+ are sinks whereas matrices in r)A3+ are sources. 
The orientation of the eigenvectors of both stability matrices are within the interval [0„iax, 7r/2]. 
Considering the latter case of spiral points, we can at least say that within one family {B/t-|? = 0...3, r — ±1} 
of matrices either Bi+ is a stable spiral and B3+ is an unstable spiral point or vice versa. 

— The class C^: Here more cases are possible. Looking at the matrices {Bj+jZ — 0...3} of the family {B;t-|Z ~ 
0...3,T = ±} there can be either 

* no spiral point, but two sinks and two sources (analogous to C^) 

* two spiral points (one stable, one unstable), one sink and one source (analogous to in C^) 

* four spiral points (two stable, two unstable) 

Which of these cases occur is a question of the variation of V'iT(0) ~ a-nd of the phase ipiTi4'nun)- If real 
eigenvalues occur, the angles of the corresponding eigenvectors are within [7r/2,0max] for matrices in ^1+ 
and A3+ ■ The position eigenvectors of stability matrices matrices in Ai+ and .43+ is not determined within 
this classification scheme. 

To determine the minimal sets Si of matrices Cka necessary for stabilisation of all FPs of a two dimensional chaotic 
dynamical system let us first consider the classes Air^iC^ whose representatives correspond to stable FPs. We herefore 
form pairs {It, Vt') abbreviating the two classes Aw and Aiit' for any C*. 



(l+,2+) or (2+,3+) 

(l+,2+) or (2+,3+) or (3+,0+) (24) 
(l+,2+) or (2+,3+) 



When looking for the minimal set of matrices Cka necessary for stabilisation one has to take into account that 
only certain kinds of FPs can occur in the system x = J^(x), eqn. (^), derived from the original dynamical system 
Xi+i = /(xi). In two dimensions /(x) can have saddle points only and the FPs of i^(x) are therefore either saddle 
points too or sinks. 

Let us discuss the set S saddle of matrices Cka which stabilise saddle points first. Since the determinant of the 
stability matrix is negative for a saddle point and positive for any stable sink or spiral point, the corresponding 
stabilising matrix Cka has the form Cfc_. Since the saddle points of all classes C^, and have to be stabilised by 
S we have to determine S such that any class Ai- of original matrices is transferred into at least one element in each 
of the pairs (l+,2+) or (2+,3+) or (3+,0+), which is the union of the pairs of all three classes C^, and in 
eqn. (p^). The transition tabic || shows that there are just two possibilities of minimal sets: Ssaddie = {Co-, C2-} 
or {Ci_, Cs-}. Each of these sets has to be combined with sets Ssink that stabilise the sinks of x = FijK.) to yield a 
possible set S. Since the latter are already stable, the identity transformation Ssink = {Co+} is sufficient. Indeed, it 
is easy to see that no other Cka is able to achieve the same: Sinks can occur in both and for the classes ^1+ 
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and A2+ and in for the classes A2+ and ^3+. We therefore can list these different classes of sinks and sets Ssink 
of matrices Cka necessary to stabilise them as follows: 



original classes 


stable classes 


sets of stabilising C 




c^r]Ai+ 




(l+,2+) 


{Co+} or {Ci+} 




ci n A2+ 


(n 


^,2+) or (2+,3+) 


{Co+}or {Ci+,C34 


-} 


n A2+ 


(n 


^,2+) or (2+,3+) 


{Co+}or {Ci+,C34 


-} 


n ^3+ 


(2H 


^,3+) or (0+,3+) 


{Co+}or {Ci+,C34 


-} 


n ^1+ 




(l+,2+) 


{Co+} or {Ci+} 




c3 n A2+ 


(n 


-,2+) or (2+,3+) 


{Co+}or{Ci+,C34 


-} 



So the smallest set that stabilises any sink independent of the classes C^, and C"^ is Ssink — {Co+}- If we now 
form the unions Ssaddie U Ssink to get a set of matrices Cka that stabilise all FPs of a two dimensional fully chaotic 
dynamical system we end up with two 'global' minimal sets: 

5i ={Co+,Co-,C2-} and ^2 = {Co+, Ci_, C3-} (26) 

There are other sets which also do the job, but they contain at least four matrices and are therefore not minimal. In 
previous numerical investigations it was observed |^,^ that the transformation belonging to the matrix Co+ yields 
a particularly large number of stabilized fixed points. From the above discussion this becomes now understandable 
since Co+ is responsible for the stabilization of the sinks in equation (|^) which occur in many different classes. 

Let us remark that the above discussion includes also sinks of the original system / since they become sinks of F 
and can therefore be treated in the same way, i.e. are 'stabilised' using the same minimal sets Si or 52- 

The concept of characterising FPs by sets of manifolds which are invariant with respect to the operations Cka 
can be extended in a natural way. One can consider manifolds which are composed by not only two subsets Zi 
and Z2 as above, but of four subsets Z1...Z4. These subsets are defined implicitly by demanding a given ratio of 
the two components of the flux. They correspond for the linear regime to four angles {^i, 02, 03, 04 } 7 which are a 
generalisation of the parameters {0min, 0max} in the discussion above. This larger set of parameters leads to a finer 
partition of the space of all 2 x 2 matrices, which in turn allows a complete assignment of stability properties to the 
different matrices. The definition of the family of a given matrix is analogous to the corresponding definition ( |l0| ) 
and reflects the action of the group of matrices {Cko-}. However, finding the proper partition is not straightforward 
and is left to a future investigation. 



III. APPLICATIONS 



In this section we apply the SD method whose theoretical background has been discussed so far. Furthermore we 
will discuss several improvements of its original algorithmic implementation to locate the UPOs |^,^ With respect 
to the detection of the orbits our aim is twofold. First we are interested in complete sets of UPOs for higher periods 
of the above maps and in analysing them with respect to the distributions of the corresponding Lyapunov exponents. 
Second we want to demonstrate the suitability of our method to detect the least unstable orbits up to high periods. 
The latter is an extension of the work given in ref. p6[ |. We remark that very recently an efficient algorithm for 
detecting UPOs in chaotic systems based on a combination of the SD method and a Newton-Raphson like approach 
has been developed and successfully apphed [p9[. 



A. Finding the Fixed Points 

Concerning the efficient algorithmic implementations of the SD method we face two main problems: 

1) The completeness of the detected set of UPOs which is of conceptual character. 

2) Separating closely neighboured UPOs which is an issue only for our particular implementation of the SD method, 
i.e. characteristic for our numerical approach. 

Let us first address the completeness problem. Of course there is no exact proof of completeness for the detected 
UPOs within the SD method. However, a properly chosen sequence of sets of initial conditions which cover the phase 
space of the dynamical system as neat as possible can significantly lower the probability of missing any UPOs. Due 
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to the presence of length scales which differ by many orders of magnitude induced by the fractal structure of the 
corresponding strange attractor we proceed here as follows. 

We introduce a set of grids Gi, i=l,2,... of starting points which are cumulative in the way that the points of Gi 
fill gaps on the attractor that are larger than a given size in the union of the preceding grids Gi L) G2 U ... U Gi_i. 
In the particular case of the Ikeda map (see below) we generate a sequence of six grids Gi,...,G6 with Gi containing 
approximately 4500 points, while Gi, i ~ 2, ...,6 contain approximately 1500 points each. The starting points of each 
set Gi are propagated with the transformed maps Sq+, So-, and 5*2-, i.e. applying the matrices Co+, Cq- and C2- 
according to (|l|). The propagation of a particular trajectory is stopped at the point x if 5 = |/^^-'(x) — x| < e, where 
e is the desired resolution of the FPs. For periodic orbits of the Ikeda map with p = lA and 15, e < 10^^° turned 
out to be sufficient to resolve the different cycles. Propagating the starting points of the grids Gi with 5o+, 6*0- and 
*S'2- yields the sets A^^ o+i -^i,o- and -/Vi,2-i * = 1, ■•■,6, respectively, of points which to good accuracy approximate 
the FPs of the Ikeda map. The set Ni = Ni^^ U Ni^g- U Ni^2- contains the FPs of the map found by propagating 
the points of Gi. Then we consider the number rii^ka of FPs that appear in a particular Ni^ka (and therefore in Ni) 
with (ka) = (0+), (0—), (2—), but which are not contained in any other Nj, j < i. If Ui^i^ka = ni_ka = 0, i.e. no 
additional FPs have been found when propagating two subsequent grids Gi-i, Gi, the set of FPs of the transformed 
system Ska is considered to be complete and the procedure of constructing subsequent grids is stopped. 

The second improvement concerns the separation of neighbouring UPOs. Using however an appropriately defined 
distance dxy between two orbits of period p, Xi = {xii,X2iY' , ifi = {yu, V2iY' , i — 1, ■ • • ,p: 

d^y^mm ^ {xj,-yj(^,+k mod p))'^ (27) 

...,p'-l l<i<p 
l<j<2 

In the following we consider a set of FPs that belong to different periodic orbits. The FPs obtained in the procedure 
of propagating the transformed systems S'o+, ^o- and S'2_ as described above provide an example for such a set. By 
forming all possible pairs of any two FPs of this set and looking at the resulting distribution of values of dxy (e.g. by 
plotting all dxy in a logarithmic scale) one can distinguish three subsets, separated by gaps differing by several orders 
of magnitude (Fig. ||). The set with the largest values of dxy contains all pairs formed with different orbits, whereas 
the second and third largest set is composed of pairs of identical orbits. A possible explanation for the appearance of 
the gap between the second and the third zone is the following: A trajectory of the transformed map approaches the 
FP on a curve which in the local but still nonlinear neighbourhood is close to the least stable of both stable manifolds 
of the FP. Two trajectories {y^} and {zj} can therefore evolve towards a FP Xq along the same or opposite directions. 
The same holds for all other FPs of the orbit f^'^\xo) with the trajectories {f^^^Yi)} and {/'^''^(z^)}, r = 1, ...,p — 1. 
In the case of an antiparallel approach to the FPs of a given cycle the contributions to dxy are much larger compared 
to the case of a parallel approach, which finally yields the gap between the second (parallel approach ) and third 
(antiparallel approach) region in Fig. ^. If the first and the second zone happen to merge, the accuracy of the FPs 
is not sufficient to distinguish different orbits and has to be refined by further propagation with the corresponding 
transformed maps. What is more, this distinction can be used to derive an estimate for the absolute accuracy of the 
FPs derived by propagating the appropriately transformed system: The maximal separation of two points belonging 
to the same cycle is given by the square root of the value of dxy at the upper edge of the lowest band in the distribution 
of values of dxy 

With the above method we selected exactly one point of each periodic orbit we found. The other points are obtained 
by simply propagating the selected point with the original map. 



B. Complete Sets of Orbits and Lyapunov Distributions 

In the following we ffist focus on the Ikeda map |^8|| . It is given by Xn+i — a + P{xn coswn — y-n sin w„), j/n+i = 
(3{xnShiWn + ynCOSWn), whcrc Wn — "f — i^x^ +y'^ ' ^^'^ attractor to be investigated appears for a = 1.0, /3 — 0.9, 
7 = 0.4 and 6 = 6.0. For the Ikeda map the periods p — 1,2, 13 have already been investigated in ref. |2^. In 
order to indicate the applicability of the previously discussed algorithmic implementation of the SD method beyond 
those periods we calculate complete sets of orbits for p = 14, 15. 

A rough hint towards the completeness of the result is the convergence of the corresponding topological entropy 



h = limp ,00 hp with hp = hin{p)/p, n(j>) being the number of FPs of all cycles of period p (see table III), h seems 

to be converged fairly well to a constant value. Our sets of orbits for the Ikeda map, together with the corresponding 
results for the Henon map, allow us to study the distributions of the Lyapunov exponents of the UPOs for the two 
different systems (see below). 
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Before entering into this discussion we comment on some observations made by applying the SD method to the 
Ikeda map. As explained in Sec. 2.2, each point of an UPO is stabilised by two matrices, which both have either 
positive or negative determinant. Therefore, we can group the points of one UPO of period p into the sets 

S'+(p) : sinks, stabilised by either Co+, Ci+, C2+, or C3+ (28) 
^-(p) : saddles, stabilised by either Cq-, Ci_, C2-, or C3_ (29) 

We found the surprising result that for all periods p = 1 — 14 both sets S'+(p) and S-{p) contain the same number of 
points, \S+{p)\ and |5_(p)|, respectively. The values are \S±{p)\ = 1,2,4,8,11,26,36,64,121,242,419,692,1262,2256 
for p = 1 - 14. Only for p = 15 we olaserve a deviation from this fact: |S'+(15)| = 4409, |5'_(15)| = 4379. At the 
current state we cannot judge whether this indicates missing orbits or represents a violation of the observed rule for 
p = 1 — 14. If we assume the validity of the equality 15+1 = \S-\ as a general law, this implies a symmetry relation 
between the orbits of the map: Each orbit with a positive determinant of its stability matrix is related to exactly 
one orbit with negative determinant. In section 3, we will seek for the most stable orbits of the Ikeda map for higher 
periods and we will find a similar pairing of orbits, i.e. two orbits, one in S+ and one in 5_, have nearly the same 
Lyapunov exponent. This might suggest a symmetry of the underlying Ikeda map. For all detected UPOs we calculate 
the Lyapunov exponent Aorb — log(|p|)/p, where p is the largest eigenvalue of the matrix M = Mp ■ . . . ■ M2 ■ Mi. 

Fig. shows the normalised distributions D{A) of Lyapunov exponents A of all orbits of order p = 24, . . . , 27 
for the Henon map |2^ given by a;„+i = 1.4 — a;^ + 0.3 yn, Un+i = Xn, and for p = 12, ... , 15 of the Ikeda map 
(primitive orbits only). For both maps the Lyapunov exponents form a band and the distributions show a more and 
more pronounced peak as p increases. The peak of the distribution appears around A = 0.5 for the Henon and around 
A = 0.68 for the Ikeda map. Around this peak, the distribution is expected generally to be well approximated by 
a Gaussian ||l2|, |Q. However, globally the distributions of both maps clearly deviate from Gaussian behaviour in 
some respects: 

• Unlike symmetric distributions, both distributions exhibit an enhancement for small values of the Lyapunov 
exponent. This indicates that there are correlations within these orbits. 

• A second feature is the occurrence of peaks in the main bulk of the distribution. A characteristic peak of Henon 
map appears at A = 0.551 and a similar spike is visible for low values at A = 0.435. 

Probably the Lyapunov distribution for the cycles of the Ikeda map displays similar features, though the number 
of cycles for period p = 15 is not enough to allow for a sufficient resolution of the distribution. 

These features show that the Lyapunov distributions contain interesting information on the systems dynamics. 



C. Stability Ordering of Cycles 

Since the SD method is, by construction, changing the stability properties of a dynamical system, it is not surprising 
that the value of the tuning parameter A is of relevance to the magnitude of the Lyapunov exponents A^JJ^^ to be 
detected. The critical parameter Xka,i of a point of a cycle stabilised with Cka is defined as the largest value of 
A for which both eigenvalues of the transformed stability matrix Ts^^ in eqn. (^) have an absolute value less than 
unity, which marks the transition from instability to stability. 

As introduced in ref. , an approximately monotonous relation between A^^J^ and the critical value A^"^ j can be 
observed. For the orbits of the Henon map, which are stabilised by a specific Cfccr~matrix, this monotonous relation 
can be clearly seen. However, when examining other maps, e.g. the Ikeda map, one finds that this relation is obeyed 
less strictly. However, a slightly different concept of ordering does the job: We consider all points rp\ i = l..p of an 
orbit j of given period p, which in general are stabilised by different CfcCT-niatrices with different critical A^"^ ^-values. 
In contrast to the approach chosen in we now allow all eight Ci:^ -matrices to be used as stabilising transformations. 
As explained in section 2.2, each point of the orbit j is stabilised by two matrices Cka, Ck'cr', with a particular A^f^^- 
To each orbit of period p there belongs a set of 2, 3 or 4 Cfco-matrices stabilising different cycle points and a set 
of 2 n values A^^lj. Now we ask for the largest Ac out of this set and call it the A^/j, of the corresponding orbit 
j. The corresponding plot for the Ikeda map is presented in Fig. |^ and shows a sufficient ordering of the stability 
coefficients of the detected UPOs with respect to the corresponding critical values X^^Ji,- The area in each subfigure 

which is shaded gray indicates the approximate position of the points A^'JJj^^x'^JJf^) for period p = 15. Obviously the 
distributions become increasingly flatter with increasing p, which corresponds to a better stability ordering of the 
respective periodic orbits. This allows us to calculate the least unstable orbits of a map in a systematic way. We 
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decrease the value of A used in eqn. (|l|) and register the stabilised points one by one. The main difference to the 
procedure in ||2^ is that we now have to consider the set of all Cfco- -matrices to find the A'^-'^ as the largest A^^2 i of all 
points Ti of an orbit. To implement these ideas, we construct a successive number of cumulative grids Gi, G2, G3, . . . 
leading to an increasingly finer covering of the attractor as described in Section 3.1. For our investigations of the 
Henon and Ikeda map we used 10 grids Gi, . . . , Gio each with about 250 points. Then we perform the following steps 
(starting with A = 0.8 and i = 1): 

1. Begin with an initial value of A and a grid Gi of points. Propagate Gi 8 times with the stabilised systems 
according to (|l|) for fixed A, using a different Cfco-, fc = 0, 3, <t = ± each time. 

2. If step 1 does not yield the desired number of orbits, replace A =^ r • A, r « 0.8 and perform step 1 again. 

3. Replace gi — > g^+i, A — > r ■ A and go to step 1. 

The procedure is converged if the set of the N most stable cycles resulting from the grid i/; = Gi U G2 U ... U Gi 
remains the same compared to the set obtained from the grid ~ Gi UG2 U ... UG^+i . As a result of this procedure 
we get a set of N orbits for each grid. 

For the two maps studied here we can clearly observe convergence. In Fig. ^ the Lyapunov exponents for the 
ten most stable orbits of each period are shown for both the Ikeda and Henon map for periods p = 1 — 36. These 
Lyapunov exponents correspond to the lower edge of the distribution of Lyapunov exponents as given in Fig. |^. 

Two features in Fig. ^ are remarkable: 

• The Lyapunov exponents of the least unstable orbits of a certain length p of both the Henon and Ikeda map are 
approximately in the same range. This is valid for all periods p considered. They both decrease with increasing 
period. This corresponds to a shift of the lower edge of the Lyapunov distribution towards decreasing values of 
A with increasing period. The inset in Fig. ^ shows the same Lyapunov exponents in a log-log plot. The mean 
of the total distribution, i.e of all UPOs for a certain period, is the average Lyapunov exponent A of the maps 
which is approximately constant. Therefore the linear decrease in Fig. || implies that the spreading W{p) of the 
tail of the distribution of the Lyapunov exponents as displayed in Fig. ^ for period p grows approximately as 
an algebraic function of the period, i.e. W{p) oc p'^ , 77 > 0. 

• Both maps show exceptionally small Lyapunov exponents ( e.g. periods p = 13, 16, 18 and p — 26, 28, 30 for the 
Henon map, periods p = 19, 21 and p = 24, 27, 30 for the Ikeda map). These orbits seem to approach the main 
part of the distribution with increasing period. 

IV. SUMMARY 

This paper has two main objectives. First it presents a novel approach towards a better understanding of the SD 
method, thereby establishing a geometric interpretation and classification, and second it provides results of applications 
of the SD method to two dimensional maps. 

In the first part we investigated the stability transformations as proposed in by introducing two classifications 
of the corresponding stability matrices of all FPs. These are based on properties which change in a regular and 
well-defined way when the stability transformations are applied. The first classification is with respect to manifolds 
which are invariant to the stability transformations. The second one mirrors the change of the azimuthal angle of the 
flux in the neighbourhood of the FPs. According to these classifications of the stability matrices we can assign to 
each set of corresponding FPs a certain dynamical behaviour related to their stability properties. This provides new 
insights into the mathematical group structure which is impressed by the stability transformations onto the stability 
matrices and FPs. Particularly this allows to determine the minimal set of transformations necessary to detect all 
FPs of a given two dimensional system. We thereby learned how simple global operations on the dynamical system 
change the stability properties of fixed points. This point of view has the advantage that it does not rely on the 
analytic expressions for eigenvalues and is therefore more suited for the extension of the method to higher dimensions. 

The second part of the present work deals with the algorithmic implementation of the SD method, its application 
to the finding of UPOs in strange attractors as well as the evaluation and interpretation of the achieved results. In 
particlar we demonstrate that even longer cycles can be detected. This is achieved by using a special sequence a 
grids of initial points for the propagation and by an improved technique to separate distinct UPOs. It becomes thus 
clear that the maximal period of orbits to be detected by the SD method is limited by the machine precision, not 
by failure of the method itself. As a result of our investigations of the Ikeda map we get distributions of Lyapunov 
exponents which show characteristic deviations from the first-order approximation by a Gaussian distribution. The 
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distributions of the Henon map show similar, but distinct characteristics. Analysing the UPOs of the Ikeda map 
for the periods p = 1 — 15 with respect to their stability in all different transformed systems suggests an underlying 
symmetry relation for this map, which implies a correlation between distinct UPOs of the same period. As a second 
numerical investigation we search for the ten periodic orbits with the smallest Lyapunov exponents for the Ikeda and 
Hcnon maps. Since these orbits are the least unstable ones, it is possible to extend the investigation up to period 
p = 36 for both maps. The distributions of the cycles as a function of the period shows a remarkably regular overall 
tendency with characteristic deviations for both maps. What is more, this part of the Lyapunov spectrum covers the 
section of small values of the Lyapunov distributions which differs most from the Gaussian approximation and might 
therefore provide valuable information on the dynamical systems. 

Finally we remark on very recent developments concerning the detection of UPOs. In ref. the SD method has 
been combined with the Newton-Raphson method in order to speed up convergence in the linear neighbourhood of 
the FP. Such hybrid algorithms are very desirable since they perform very efficiently while preserving the desired 
global character of the SD method. 
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Table Captions 



TABLE I. Multiplication table for the matrices Ctcr 



TABLE IL Transition matrices Cka necessary for transitions between different classes Aw of stability matrices of fixed 
points. The first three columns indicate the combination of classes occurring in a two dimensional chaotic dynamical system. 



TABLE III. The number of prime cycles with period p, the total number of cycle points of order p and the topological 
entropy hp for Ikeda map, periods p = 12, 15 



Figure Captions 



FIG. 1. Definition of the polar angles cf> and 'ipii^i{(j>) of the displacement x = (cos 0, sin (^)^ relative to the FP at (0,0) and 
the flux Ffeo-(x) — (cos , sin i/);/^/)"^. (ka) indicate the particular stability transformation applied according to eqn. 
The indices {I't') are given by eqn. (^) using (ka) and the indices (It) of the class Air to which the stability matrix of the 
original FP belongs. 

FIG. 2. a) Phase portraits of fixed points with det(B;T-) > 0, < 4>min, 4>max < tt/2. b) and c) are the corresponding 
i/'!T(</')-diagrams. In a) the manifolds Zi, Z2 (long dashed lines), eigenvectors (lines with filled arrows) and some trajectories 
are shown. The areas shaded gray indicate the intervals of the locations of the eigenvectors. They correspond to the shaded 
boxes in the diagram c) which show these intervals for the four FPs discussed. 



FIG. 3. same as Fig. 2, but with det(BiT) < 0, < (l>min, (f)max < tv/2, 



FIG. 4. same as Fig. 2, but with det(B;T) > 0, < (f'min < vr/2 < 0, 



FIG. 5. same as Fig. 2, but with det(B;T) < 0, < 4>min < 7r/2 < (l>max, 

FIG. 6. Distribution of the distances d^y of the orbits of the Ikeda map for period p = 15 of the set A^i. The abscissa of 
the figure covers the 16076 pairs of FPs in A^i. 

FIG. 7. Normalised distributions of the Lyapunov exponents (prime cycles only) of the Ikeda and Henon maps for various 
periods. 

FIG. 8. Distribution of the Lyapunov exponents of the stabilised orbits as a function of the critical stabilising parameter 
\i'X for Ikeda map, periods p = 10 — 15. The areas shaded gray show the approximate range of the distribution for p = 15. 

FIG. 9. The Lyapunov exponents of the 10 least unstable orbits of Ikeda and Henon map for p — 1...36 The inset shows the 
same distribution on a log-log scale. 
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Table I 
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Table III 



period p 


12 13 14 15 


no. of prime cycles 


110 194 317 576 


no. of fixed points 


1384 2523 4512 8668 


topol. entropy 


0.603 0.603 0.601 0.604 



